% 使用conv计算小波边缘提取
% P 要提取的图像
% scale 最大尺度
% threshold 阈值，常数，需要修改
% wavelets_index 使用哪种小波
% E 返回的边缘提取后的边缘图像
% M 返回的边缘提取后的模值图像
function [E M] = compute_wavelets_conv( P, scale, threshold, wavelets_index )

% 说明:
% 1 小波变换使用 a'trous 快速小波提升算法(即教材P96 P188提及的 多孔算法) 见参考文献 《二维A Tuous算法图像边缘检测研究》;
% 2 二维二进小波变换快速算法的推导过程 见教材P195  式(7-42)是推导结论;

% 目前已解决问题:
% 1 二维二进小波滤波器的计算构造(实质是求两尺度方程系数 目前没有成熟的计算方法 
%   仅能查到结果并且有纰漏 选取教材提供的两组数据作为解决方案)
%     1.1  非正交二次样条小波
%     1.2  正交二次样条小波(近似解与修正解)
% 2 在确定阈值T的情况下边缘检测
%     2.1  根据滤波器系数进行小波变换求像素点的梯度值和梯度角
%     2.2  非极大值抑制
%     2.3  阈值化
%     2.4  边缘连接

% 目前遗留问题:
% 3 阈值的自适应算法 或者通过GUI作为手工输入参数处理
% 4 多个尺度边缘下结果的筛选处理
% 5 工程问题: 图形界面中应该设置哪些参数及实现;

I=P;
[N,M] = size(I); 

%以下h g delta三个变量的意义: 分别对应于教材P197(7-42) 的 h g l滤波器系数

%1 非正交二次样条小波系数的构造 参考 教材P182例7.1和P193例7.4 需要阈值化得到边缘
% h = [0.125,0.375,0.375,0.125]; 
% g = [-1.0,1.0]; 
% delta = [0.015625,0.09375,0.234375,1.3125,0.234375,0.09375,0.015625];

%2 正交二次样条小波系数的构造 参考 教材P183例7.2和P194例7.5 需要阈值化得到边缘
%h = [0.125,0.375,0.375,0.125]; 
%%g = [0.0061,0.0869,0.5798,-0.5798,-0.0869,-0.0061];%近似解 
%g = [-0.9324,0.6736,0.2253,0.0335];%修正解 from PPT

if wavelets_index==1
	% 正交二次样条
	h = [0.125,0.375,0.375,0.125]; 
	g = [-0.9324,0.6736,0.2253,0.0335];%修正解 from PPT
	delta = [0.00003,0.00727,0.03118,0.06623,0.79113,0.06623,0.03118,0.00727,0.00003];
elseif wavelets_index==2
	% 正交三次样条
	h = [0, 0, 0.0625, 0.2500, 0.3750, 0.2500, 0.0625, 0, 0];
	g = [0.00008, 0.01643, 0.10872, 0.59261, 0, -0.59261, -0.10872, -0.01643, -0.00008];
	delta= [ 0.0035 0.0276 0.0938 0.1795 1.1098];
else
	% 非正交二次样条
	h = [0.125,0.375,0.375,0.125]; 
	g = [-1.0,1.0]; 
	delta = [0.015625,0.09375,0.234375,1.3125,0.234375,0.09375,0.015625];
end

%参数K意义:多尺度边缘检测中的尺度个数
K = scale;
%tan_angle = zeros(N, M); %角度矩阵(此处采用角度，也可直接用正切来比较)

%四维矩阵  实际上只使用到其中的后两维存储二维小波变换的结果
a(1:N,1:M,1,1:K+1) = 0; 
dx(1:N,1:M,1,1:K+1) = 0; 
dy(1:N,1:M,1,1:K+1) = 0; 
d(1:N,1:M,1,1:K+1) = 0; 
tan_angle(1:N,1:M,1,1:K+1) = 0;


%边缘化第一步 对I进行处理
a(:,:,1,1) = conv2(h,h,I,'same'); 
dx(:,:,1,1) = conv2(delta,g,I,'same'); 
dy(:,:,1,1) = conv2(g,delta,I,'same'); 

x = dx(:,:,1,1); 
y = dy(:,:,1,1); 

d(:,:,1,1) = sqrt(x.^2+y.^2); %计算1尺度下的梯度模值
dt=d(:,:,1,1);%dt at为中间变量
dt2=d(:,:,1,1); % 存放没有进行极大值提取的边界图像
tan_angle(:,:,1,1) = atan(x./y);%计算1尺度下的梯度角
at=tan_angle(:,:,1,1);


%沿着梯度方向进行 非极大值抑制
%dt在非极大值抑制的过程中会发生改变故将原始dt存于PT 然后基于PT对dt进行非极大值抑制
dt = do_modmax( dt, at );

%阈值化
P=zeros(N,M);
for i=1:N
	for j=1:M
		if dt(i,j)>threshold
			P(i,j)=255;
		end
	end
end

%边缘连接
P = do_edgeconnect( P, at );

lh = length(h); 
lg = length(g); 

%接着进行J个尺度下的边缘检测 每次对上一次的低频a(:,:,1,j)进行处理
for j = 1:K
	lhj = 2^j*(lh-1)+1; 
	lgj = 2^j*(lg-1)+1; 
	hj(1:lhj)=0; 
	gj(1:lgj)=0; 
	for n = 1:lh 
		hj(2^j*(n-1)+1)=h(n); 
	end

	for n = 1:lg
		gj(2^j*(n-1)+1)=g(n); 
	end

	a(:,:,1,j+1) = conv2(hj,hj,a(:,:,1,j),'same');
	dx(:,:,1,j+1) = conv2(delta,gj,a(:,:,1,j),'same'); 
	dy(:,:,1,j+1) = conv2(gj,delta,a(:,:,1,j),'same'); 

	x = dx(:,:,1,j+1); 
	y = dy(:,:,1,j+1); 

	% 计算j+1尺度下的梯度模值
	dj(:,:,1,j+1) = sqrt(x.^2+y.^2); 
	dt=dj(:,:,1,j+1);
	dt2=dj(:,:,1,j+1);
	% 计算j+1尺度下的梯度角
	tan_angle(:,:,1,j+1) = atan(x./y);
	at=tan_angle(:,:,1,j+1);

	%沿着梯度方向进行非极大值抑制
	dt = do_modmax( dt, at );
	%阈值化
	P=zeros(N,M);
	for i=1:N
		for j=1:M
			if dt(i,j)>threshold
				P(i,j)=255;
			end
		end
	end
	%边缘连接
	P = do_edgeconnect( P, at );
end

E = uint8(P);

Emax = max(max(dt2));
M = uint8(dt2*255/Emax); %TODO

end

